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Abstract. The effects of a single non-magnetic impurity on superconducting states in the 
Penson-Kolb-Hubbard model have been analyzed. The investigations have been performed within 
two steps: (i) the homogeneous system is analysed using the Bogoliubov transformation, whereas 
(ii) the inhomogeneous system is investigated self-consistent Bogoliubova-de Gennes equations 
(by exact diagonalization and the kernel polynomial method). We analysed both signs of pair 
hopping, which correspond to s-wave and r]-wave superconductivity. Our results show that an 
enhancement of the local superconducting gap at the impurity-site occurs for both cases. We 
obtained that Cooper pairs are scattered (at the impurity site) into the states which are from the 
neighborhoods of the states, which are commensurate ones with the crystal lattice. Additionally, 
in the 77 -phase there are peaks in the local-energy gap (in momentum space), which are connected 
with long-range oscillations in the spatial distribution of the energy gap, superconducting order 
parameter as well as effective pairing potential. Our results can be contrasted with the experiment 
and predicts how to experimentally differentiate these two different symmetries of superconducting 
order parameter by scanning tunneling microscopy technique. 
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1. Introduction 

The Penson-Kolb-Hubbard (PKH) model is one of the conceptually simplest phenomenological 
models for studying correlations and for description of superconductivity (SC) in narrow-band 
systems with short-range, almost unretarded pairing. Generally, impurities destroy or worsen 
desirable SC properties, but sometimes can induce new interesting phenomena. 

The cuprates and iron-pnictides are examples of superconductors with short coherence length. 
Moreover, such systems are highly inhomogeneous and the role of impurities on SC could be crucial. 
In both groups of materials the superconducting state itself is generated by chemical doping which 
inevitably disorders the samples, and second, local probes of the quasi-particle states near the 
impurity sites can provide important information on the underlying system [1, 2]. In discussed 
materials the physical properties changing by doping, the result of which are phase transition. 
The phase diagrams for these materials include antiferromagnetically ordered (AF), metallic (non- 
ordered) and superconducting (SC) phases. With increasing concentration of current carriers the 
AF phase is vanishing first, then the SC phase occurs in define range of doping [3, 4]. The 
other interesting feature of these compounds are effects introduced by (nonmagnetic) impurities 
such as e.g. inhomogeneity superconducting state [5, 6, 7], pining of vortexes [7, 8], modification 
persistent current in superconducting ring [9], and induced spin density waves [10]. Notice that the 
effects of nonmagnetic impurities in above materials are qualitatively different than in conventional 
superconductors, any perturbation that does not lift the Kramers degeneracy of these states does 
not affect the mean-field superconducting transition temperature. 

The pair-hopping term (</) was proposed in Ref. [11, 12] and can be derived from a general 
microscopic tight-binding Hamiltonian, where the Coulomb repulsion may lead to the pair hopping 
interaction J = (ii\e 2 / r\j j) [13, 14, 15]. In such a case J is positive ( repulsive model J > 0, favoring 
rj-wave SC), but in this case the magnitude of J is very small. However, the effective attractive form 
(J < 0, favoring s-wave SC) is also possible (as well as an enhancement of the magnitude of J > 0) 
and it can originate from the coupling of electrons with intersite (intermolecular) vibrations via 
modulation of the hopping integral or from the on-site hybridization term in the general periodic 
Anderson model (cf. e.g. [4, 16, 17, 18, 19] and references therein). It can also be included in the 
effective models for Fermi gas in an optical lattice in the strong interaction limit [20]. The role 
of J interaction in a multiorbital model is of a particular interest because of its presence in the 
iron pnictides [21]. Notice that it has been found that SC originating from the J > 0 interaction is 
unique in that it is robust against the orbital (diamagnetic) pair breaking mechanism [19]. 

In this paper we investigate a role of single nonmagnetic impurity on SC states ( s-wave as well 
as i]-wave). We obtain the spatial distribution of local energy gap (and its Fourier transform), 
superconducting order parameter as well as effective pairing potential. We also propose an 
experimental method that can determine what type of superconductivity ( s-wave or r/-wave) occurs 
in the material. 

2. Model and technique 

We explore two dimensional square lattice of a size N x x N y with the periodic boundary conditions, 
where we assume the possibility of the hopping of (i) electrons (with amplitude t) and (ii) pairs of 
electrons (with amplitude J) between nearest neighbors (NN). Electrons with opposite spin at the 
same site interacts with energy U (inter-site Coulomb interaction). In the model we also introduce 
the diagonal disorder induced by non-magnetic impurities. The Hamiltonian of the system in real 
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space can be describe as H = H 0 + H int , where: 

H 0 = + (!) 

Hint = U Y1 c it c l| c 4 c it + J Y1 4t c 4 c 4 c it- (2) 

i (».j> 

Yl(ij) restricts the summation to NN, Cj CT (c| CT ) are annihilation (creation) operators of electron at 
*-th site with spin a = 4 is chemical potential and V) is non-magnetic impurity potential at 

*-th site. The electron hopping amplitude (t) will be taken as a scale of energy in the system. 



Figure 1. Schematic illustration of the model analyzed. The impurity can be identified as a 
change of a local potential at the impurity site (not shown explicitly). 

In the mean-field Hartree-Fock approximation [4] interaction Hamiltonian (2) takes a form: 

H“( = UJ2 (A*c 4 c it + H.c .) - |A 4 | 2 (3) 

i i 

+ J J2 ( A i c 4Cit + H.c.) a;.A,- 

(i<o) (i,j) 

where we introduce dimensionless superconducting order parameters (SOP) Aj = (c*Thus, in 
the non-homogeneous at every site the effective (local) energy gap exists, which value is determined 

by 

v i = u± i + jY,*i- (0 

(ij) 

This quantity is connected with the effective pairing potential defined as UfH = r D i /\ Aj|, which 
is also not a constant value in a heterogeneous system. To analyze the gap distribution in the 
momentum space we introduce also its Fourier transform (FT) defined as 

Vk = 

x y i 

k denote a momentum of Cooper pair. 

The further analyses are performed in two steps: (i) the homogeneous system is analyzed 
using the Bogoliubov transformation, whereas (ii) the inhomogeneous system is investigated by self- 
consistent Bogoliubova-de Gennes equations in real space approach [22], (by exact diagonalization 
and the kernel polynomial method [23]). The both approaches will be described briefly below. 
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2.1. Homogeneous system 

In a general case the SC state is characterized by the formation of the Cooper pairs with total 
momentum Q. In a case Q / 0 the SC state will be called Fulde-Ferrell-Larkin-Ovchinnikov 
(FFLO) state [24, 25, 26]. Notice that Q = 0 = (0,0) case (favouring by J < 0) corresponds to 
s-wave SC (BCS-like SC), whereas Q = II = (n, 7 r) (favouring by J > 0) corresponds to rj-wave 
SC [16], in which the SOP alternates from one site to the neighbouring one. In particular, for square 
lattice one can distinguish two sublattices in which the phase of SOP will differ by 7 r, whereas e.g. 
for triangular lattice one can distinguish three sublattices in which the phase of SOP will differ by 
27 t/ 3 [27]. In the homogeneous system only the component V k (Eq. (5)) with k = 0 or k = II is 
nonzero (s- or rj-wave , respectively). Thus, in homogeneous system (V) = 0) the SOP can be derived 
as Aj = Ao exp(iQ ■ rf). Moreover, it has been shown that flux quantization and the Meissner effect 
appear in s-wave as well as in rj-wave SC state [28]. The mean-field Hamiltonian H hf = + 

in the momentum space takes the form [29]: 

h hf = {_ tlk _ ^ c { a c ka ( 6 ) 

key 

+ U e Q ff ]T (AgC_fc + Qj,Cfc t + H.c.) - Uq^N\A 0 \ 2 

k 

where 7 k = 2 (cos(k x ) + cos(k y )), U k * = U + Jy*, is a effective pairing interaction in the 
momentum space, and N = N x x N y . Moreover, in this case . Using the Bogoliubov 

transformation [30, 31] one can obtain a spectrum of the Hamiltonian as 

£fc,± = (E k | - E_k+Ql) /2 (7) 

± ^(E kt + E- k+Ql ) 2 /4+\U e Q ff A 0 \2 

Straightforward calculations lead to the following form of the grand canonical potential S2 = 
—kgT In Tr[exp(—/3H hf )]\ 

U = - k B T ^ ln(l + exp(-/3£ fcT )) ( 8 ) 

fc,T=± 

- U e ^N\A 0 \ 2 + Y,E k ,- 

k 

where /? = 1 /fcgT. Notice that Eq. ( 8 ) is also valid for the inhomogeneous system, but in that 
case the Hamiltonian spectrum can be determined only numerically. The ground state for fixed 
model parameters is founded by a minimization of f l with respect to Ao- It enables, among other, 
a founding the ground state phase diagram (cf. Fig. 2), from which we choose values of model 
parameters to further numerical calculation in inhomogeneous system. 


2.2. Inhomogeneous system 


The superconducting states in an inhomogeneous system are described by the following mean-field 
single particle Hamiltonian written in Nambu space: 


H = 


(i,j) 


4 t Ci i 


)( 


Hi] t 
A*. 

i] 


- h: 


ill 


c it 

r t 

c ii 


( 9 ) 
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Figure 2. The ground state phase diagram of the model in homogeneous system (V) = O.Ot and 
// = O.Ot). White dashed lines show (indicatively) boundaries of the phases. 


The off-diagonal elements A* j = <5jj(Z7Aj + J ^ A*,) are the on-site SOP at the site i, while the 
diagonal elements f/, 7(T = —tij + (Vi — /j,)) 5ij are the single particle Hamiltonian in the presence of 
a given configuration of impurities. 

The inhomogeneous system described by Eq. (9) can be solve using the Bogoliubov-Valatin 
transformation: 


Cic t — 


£(* 


( 10 ) 


where 7 na and r 'fj la are the quasi-particle operators, Ukna and Vkna are the Bogoliubov-de Gennes 
(BdG) eigenvectors. One can obtain the self-consistent BdG equations in real space: 


= £ 


Hi. 


A ij 
- H* 


*jncr 


(ii) 


where and A,; 7 have been defined previously. The on-site SOP are given by: 


A i — (c, ; C ; y ) 

= 'y "j ( u in^ v inlf (^n^) _ u inl v in^f ^ni)) > 
n 


( 12 ) 


where f(E) = 1/ (exp(/3E) + 1) is the Fermi-Dirac distribution function. 

This method enables to determine the SOP in inhomogeneous system in self-consistent 
way [5, 6, 9, 32, 33, 34, 35, 36, 37, 38]. Because it is necessary to solve eigenproblem (11) of 
the Hamiltonian, the method can be used for systems of a size of the order 40 x 40. Thus we 
also implement kernel polynomial method [23], which is the extension of Bogoliubov-de Gennes 
equations in Chebyshev polynomial basis [39, 40, 41, 42, 43, 44, 45, 46]. Such a method is iterative 
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and enables to determined SOP in systems, which size is limited by computing time. Below we 
briefly discuss the main ideas of this approach. 

The Chebyslrev polynomials can be written as < j> k (x) = cosffc arccosz] and satisfy recurrence 
relation: 


4>o{x) = 1, 4>l(x) = X, 

4>k+i(x) = 2 x<j> k {x) - <j) k - 1 (®), (k > 1). 

The polynomials <f> k (fc = 0,1,...) are orthonormal basis with interval [—1,1]: 

6{x - x') = y ^ k (x)<j>k(x'), x € [-1,1] 

Wk 


W{ x) = 


n —0 

l 


\J\ — X 1 


W k = 


tt( 1 - 4o) 


(13) 


(14) 


(15) 


To obtain eigenvectors of Hamiltonian (9) in the basis of function (j> k one need to scale its matrix 
representation: 

H' = (H — bl)/a, e 7 = (£ 7 — b)/a (16) 

^ {£max £min) / 2, b (£max T £min) / 2, 

where £ min < 4 < £ 


min _ max 


is the range of eigenenergies of hamiltonian H 1 H' is a scaled Hamiltonian 
with (also scaled) dimensionless eigenvalues e 7 g [—1,1]. Using the constant 2A( c A^ / -dimensional 
real vectors [e(i)] 7 = Si i7 and [h(i)] 7 = <5j+jv a .jv v , 71 we obtain 


A i = (cqc it ) = e T (i)h n (i) 


n —0 


T 




where 


T n = J dxf(ax+ b)W(x)(/) n (x ), 

e„(*) = (fn(H')e(i ), ft.„(i) = cf n (H')h(i). 

A sequence of the vector = cj) n (H')q (q = e(i), h(i)) is recursively generated by: 
Qn+ 1 = 277 q n Qfn — 1 (^7 'll 2), 

<?i = M H ')q, q 0 (H') = q. 


(17) 


(18) 


(19) 


T n does not depend on the index *. Therefore, the calculation of T„ can be done before any other 
calculations. The use of the recurrence formula leads to a self-consistent calculation of the BdG 
equations, without any diagonalization of 77. 

The sum in (17) is truncated to a certain cutoff M and such an abrupt truncation generally 
results in unwanted Gibbs oscillations. For this reason all coefficients T n must be multiplied by 
damping factors g n . A possible choice of the damping factors can be found in Ref. [23]. 


3. Results and discussion 

We analyze the two dimensional square lattice with the periodic boundary conditions at kT = 10 -5 f 
and p = O.Ot (it corresponds to electron concentration n = 1 independently of the other model 
parameters). 
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3.1. Homogeneous system 

Using the mean-field approach described before, we obtain the J/t vs. U/t phase diagram for 
homogeneous system (V) = 0) of a size 600 x 600 shown in Fig. 2. The diagram is nonsymmetric 
with respect to J = 0 and consists of three regions, separated by dashed lines, in which different 
phases occur. The p-wave SC can occur only for J > 0, whereas s-wave phase can be stable for 
J < 0 as well as for J > 0 (in restricted range). A necessary condition for SC phases occurrence 
is < 0 (or < 0), thus the regions of SC phases must be restricted at least by lines 

U ± 4J = 0. However, the SC (s- or rj-wave) phase is stable only if \Uf^\ is higher than some 
critical value and the boundary of stability of the particular phases determined by minimization of 
fl, denoted by dashed line on the diagram, are moved towards these lines such as shown in Fig. 2. 
On the diagram the value of SOP |Ao| is also shown. It should be stressed that computations in 
the real space and the momentum spaces give the same results. 

3.2. Inhomogeneous system 

Next we analyze the system of the size 40 x 40 with a single impurity Vi = Vi mp = 2.0 1 located 
at the center site of the system with coordinates (N x /2, N y /2) (Vi = 0 at other sites). In Figs. 3 
there are shown the dependencies (as a function of lattice site) of local effective energy gap "D, 
(Eq. 4). The panels in these figures, respectively for s-wave and rj-wave SC, are obtained for the 
same (approximately) distance from the boundaries with non-ordered phase shown (for a given 
phase) in Fig. 2 (denoted as dashed lines). It corresponds to constant value of U e ff at a particular 
phase. 

The \T>imp\ (|T>j| at the impurity site) initially for attractive U is smaller than \Doo\ (\T>i\ far 
from the impurity). It increases with increasing U and finally \Di m p\ > |2?oo| for repulsive U. The 
valley at the impurity site vanishes faster in s-wave phase, whereas in the rj-wave phase it extends 
over a larger space (over larger number of sites). Notice also that spatial variability of I'D, is bigger 
in g-wave phase than in s-wave phase (cf. also Fig. 4 for the system size of 150 x 150). These effects 
are better visible for the larger sizes of the system and is discussed in details below. 

From Eq. (5) one can obtain the distribution of the effective energy gap in the momentum 
space. In the homogeneous system with the SOP of s or r/ type all Cooper pairs have the same 
momentum Q = 0 or Q = II, respectively, and only with k = Q are nonzero. In the presence 
of the impurity V , /mp , the Cooper pair are scattered at the impurity site and thus in the system 
there are also pairs with Q ^ 0,11 (cf. Fig. 4 panels c and d). In particular, for the rj-wave 
case, there are pairs with momentum near Q = 0 without distinguished direction (central peak at 
Fig. 4.d). Moreover, one can also separate pairs with momentum in [1, ±l]-direction (four peaks) 
located near the central peak. These peaks are associated with long-range oscillations of gap in real 
space [1, ±l]-direction (Fig. 4.b). The peaks of at [±7r,±7r] still exist, but they are not shown 
in Fig. 4.d explicitly, because their magnitudes are much bigger. In the s-wave phase we observe 
the scattering at the impurity site of the Cooper pairs with Q = 0 into states with momentum near 
Q ^ 0 without distinguished direction (Fig. 4.c). Moreover, in the presence of the impurity in the 
system, 2?*. has also four peaks at [±7r, ±7r], what means that the electrons associated with rj-wave 
superconductivity become an important component of pair states. Notice, the peak of T>k at [0, 0] 
still exists, but it is not shown in Fig. 4.c explicitly, because its magnitude is much bigger. 

It is also interesting to discuss the spatial distribution of the SOP A; and effective pairing 
potential U e ff in both cases, although they can not be experimentally measured in contrary to the 
local effective gap. The SOP is reduced at the impurity site in the both cases (Figs. 4.e-f). The 
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behavior of f7 e // is rather unusual, because its value is enhanced at the impurity site (Figs. 4.g-h). 
In is also worth nothing that similar long-range oscillations of these both quantities (i.e. A* and 
U e ff in the [1, ±l]-direction in the real space (similarly as these for Vi) are far more visible in the 
p-phase than in the s-wave phase. 

The different values of the on-site impurity potential Vi mp has been also investigated. However, 
the value of Vi mp has no influence on the main qualitative results of the paper. The long-range 
oscillations of the energy gap (an other quantities analyzed) in real space (in the [ 1 , ± 1 ]-direction) 
occurs in the rj-wave phase, whereas in the s-wave phase they are absent. Only the quantitative 
changes are introduced by varying the value of V) mp (the changes of magnitudes of the oscillations). 

Contrasting our results with the effects of impurities on superconducting the FFLO phase 
has interesting results. In this phase Cooper pairs have non-zero total momentum. It is generally 
believed that that phase is very sensitive to inhomogeneities [47] because of the scattering of Cooper 
pairs on impurities. However, some kind of disorder in the system can lead to stabilization of FFLO 
phase [32], The off-diagonal disorder can also lead to local increase of SOP at impurity. 

4. Summary and final remarks 

In this report we have analyze the influence of the nonmagnetic impurity on superconducting 
properties (both s- and 77 -wave) in the materials with local electron pairing. We have shown that 
the scattering of Cooper pairs is into to the states which are from the neighborhoods of the states 
corresponding to the orderings commensurate with the crystal lattice. Additionally, in the 77 -phase 
there are peaks in the (Fourier transform of) local superconducting gap, which are connected with 
long-range oscillations of the local energy gap, superconducting order parameter as well as effective 
pairing potential in the real-space distribution of these quantities (cf. also with Ref. [48]). It is also 
interesting that the energy gap is enhanced at the impurity site for sufficiently large J, even if the 
on-site interaction U is repulsive. 

Notice that we do not consider the charge and magnetic orderings. They can occur for U > 0 
[16]. The more realistic analyses of the impurity impact on these states will be a topic of further 
works. 

The results presented in this paper can be verified, for example, by scanning tunneling 
microscopy (STM) technique. This technique can be used to study impurity states in 
superconductors. As a first test of theories, this allows a direct comparison of local electronic 
features in tunneling characteristics with the theoretical predictions for the density of states and 
superconducting local gap (for review see e.g. Ref. [1, 6 ] and references therein). Thus, in 
correspondence with the our theoretical results for the local effective gap, the STM spectroscopy 
can be useful to distinguish s-wave and rj-wave superconductivity in real materials. The results of 
this paper can give a proposal how to differentiate these two phases experimentally. 
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a. |Z>.| for U=-2.0t, J=—0.625f 



C. |Z>,| for E7=-1.0f , J=—0.875< 



e. \D t \ for U =0.0f , J=-1.125f 



g. |X>,| for I7=1.0f, J=—1.375f 
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Figure 3. Energy gap [Di\ for chosen U and J in the s-wave and y-wave phases (left and right 
panels, respectively) for n = O.Ot and Vim V = 2.0 1 (the system of a size 40 X 40). 
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Figure 4. Plots of energy gap \T>i\ (a and b), its Fourier tranform I'D*.I ( c and d), 
superconducting order parameter A* (e and f), and effective pairing interaction U (g and 
h) for s-wave (U = 2.0 1 and J = —2.0 1) and rj-wave phase (U = 2.0 1 and J = 2.0 1) - left and 
right panels, respectively. Results for the system of a size 150 X 150. 






















